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Research  in  Climatology 

Foreword 

The  Courant  Institute  of'  Mathematical  Sciences  began  a  major 
research  program  in  climatology  under  ARPA  sponsorship  through  Grant 
DA-AR0-D-31-124-72-G113  for  the  period  Feb.  1^,    1972  -  May  31,  I973 
and  the  work  continued  under  Grant  DA-ARO-D-31-124-73-G150  for  the 
period  June  1,  1973  -  May  3I,  197'^.   The  principal  objectives  were 
to: 

a)  determine  the  reliability  of  climatic  predictions  that  result 
from  the  long-term  numerical  solution  of  general  atmospheric 
circulation  models  (GCM), 

b)  devise  more  efficient  numerical  methods  for  solving  the  initial 
value  problem  for  partial  differential  equations  that  are 
defined  over  a  rotating  sphere, 

c)  develop  a  theoretical  basis  for  analyzing  the  behavior  of 
solutions,  for  large  time,  of  systems  of  equations  of  a  type 
that  model  atmospheric  flows,  and 

d)  examine  the  solutions  of  zonally  averaged  atmospheric  models. 

The  achievements  of  the  research  team  are  described  in  the 
final  report  on  the  above  grants,  by  its  two  principal  investiga- 
tors. Prof.  J.  J.  Stoker  and  Prof.  Eugene  Isaacson. 

The  ARPA  sponsorship  continues  under  Grant  AFOSR- 7^-2728  to 
extend  the  above  program.  ■■ 


P.  D.  Lax 
Director 


1.  Introduction 

The  work  in  climatology  at  the  Courant  Institute  of 
Mathematical  Sciences  of  New  York  University  (CIMS-NYU)  began  in 
February  1972  under  ARPA  sponsorship.   This  is  a  final  report  on 
the  activities  credited  to  Grant  DA-AR0-D-31-124-72-G113  for  the 
period  Feb.  I5,  1972  -  May  ^1,    I973  and  to  Grant  DA- ARO-D-5 1-124- 
73-GI50  for  the  period  June  1,  1973  -  May  3I,  1974. 

Of  primary  interest  is  the  determination  of  the  reliability 
and  efficiency  of  schemes  for  numerical  integration  over  a  long 
time  period,  of  a  system  of  partial  differential  equations  for  fluid 
motion  on  a  sphere.   That  is,  the  work  is  intended  to  contribute  to 
the  understanding  of  the  climatological  significance  of  a  long  term 
numerical  integration  of  an  atmospheric  global  circulation  model 
(GCMj  such  as  that  of  Mintz-Arakawa  [7  ;10;4  ].    The  simplest 
mathematical  system  for  fluid  flow  over  a  sphere  is  that  of  a  one 
layer  incompressible  model.   In  Section  2,  the  equations  are 
presented.   They  have  been  extensively  used  to  study  the  efficiency 
and  accuracy  of  numerical  integration  methods  for  forecasting 
purposes  [3  ;5  ;6  ;13].   The  unknown  variables  are  essentially  the 
height,  h,  of  the  layer  and  the  two  tangential  velocity  components, 
u  and  V.   A  tangential  friction  force,  proportional  to  the  square 
of  the  velocity,  is  introduced  into  the  two  velocity  equations,  in 
order  to  represent  the  dissipation  as  in  the  Mintz-Arakawa  model 
(see  [4]). 


Admittedly,  a  GCM  would  have  to  treat  the  oceanic  circulation,  if 
the  model  were  to  be  used  for  truly  long  period  numerical  studies. 
The  GCM's  of  Mintz-Arakawa,  of  NCAR,  of  NOAA,  and  of  NASA-GISS  do 
not  at  present  incorporate  a  model  for  the  oceanic  circulation. 


In  order  to  simulate  the  climate,  it  -was  proposed  to  consider 
a  flow  that  could  be  described  approximately  as  follows:   There  is 
a  low  pressure  center  at  each  pole  and  three  "semi-permanent"  high 
pressure  cells  in  each  hemisphere,  uniformly  and  symmetrically 
placed  at  about  ±30   latitude.   Denote  this  basic  climatological 
flow  by  (u  , V  , h  )  =  W  .   The  differential  equations  for  ¥  =    (u, v,h) 
are  represented  by 

(1.1)  L(¥)  =  L(W^)  , 

with  W(0)  =  ¥  (O)  as  the  initial  conditions. 

Hence,  a  solution  is  given  by  ¥  =  ¥  .   But,  the  primary 
objective  of  this  work  is  to  investigate  the  numerical  solution  of 
(l.l);   that  is,  to  find  the  solution  ¥.  of  a  finite  difference 
approximation  of  (1.1): 

(1.2)  ^A^^A^  "  -^^^c^  '■   ^A^°^  =  ^c^°^  • 

In  (1.2),  L  represents  a  discrete  approximation  of  the  differential 
operator  L. 

Three  different  numerical  methods  (i.e.,  L  )  have  been  pro- 
grammed for  the  approximate  solution  of  the  system  of  equations 
(1.1).   The  friction  coefficient  is  chosen  so  that  the  "spin-down" 
time  of  the  homogeneous  form  of  system  (1.2)  is  about  10  days.   The 
schemes  are  described  in  Sections  3,  4  and  5. 

In  order  to  represent  the  basic  climatological  flow,  ¥  ,  a 
study  was  made  of  the  motion  of  vortices  in  a  planar,  one  layer 
model.   H.  J.  Stewart's  work  [11]  on  the  stability  of  certain  planar 


vortex  "  solutions"  was  extended  by  G.  Morikawa  and  E.  Swenson 
[  9 ;  8  ] .   That  is,  Stewart  used  three  anti-cyclonic  vortices  situa- 
ted on  a  circle  of  radius  r,    to  model  the  three  semi-permanent, 
sub-tropical  highs  of  the  northern  hemisphere,  while  Morikawa  and 
Swenson,  in  addition,  introduced  a  cyclonic  vortex  at  the  center  of 
the  circle  to  simulate  the  effect  of  the  polar  low.   In  preparation 
is  a  report  by  L.  Bauer  and  G.  Morikawa  on  the  choice  of  parameters 
that  yield  climatologically  realistic  motions  of  the  vortices, 
along  with  a  study  of  the  stability  of  the  system  for  the  complete 
range  of  parameters.   The  vortex  approximation  of  the  stationary 
highs  and  the  polar  lows  is  singular.   The  mathematical  singulari- 
ties are  removed  smoothly  to  generate  the  basic  flow,  ¥  .   A  brief 

account  of  this  vortex  work  leading  to  the  definition  of  W  is 

°  c 

given  in  Section  2. 

In  order  to  provide  a  more  realistic  approximation  for  the 
basic  spherical  flow,  the  study  of  spherical  vortex  flows  was  begun 
by  S.  Friedlander  and  A.  S.  Peters.   S.  Drew  analyzed  the  general 
circulation  pattern  that  results  from  a  zonally  averaged  model. 
M.  Ghil  completed  a  preliminary  study  of  the  initialization  problem. 
A  listing  is  given  in  Section  6  of  all  of  the  papers  and  reports 
credited  to  this  program,  along  with  a  brief  account  of  the  ongoing 
work  that  has  not  yet  been  published. 

In  Final  Report,  II  for  the  ARPA  Grant  AFOSR- 7^-2728,  exten- 
sions of  this  climatology  activity  will  be  described. 


2.  The  Mathematical  Model,  Coefficient  of  Friction,  Zonal  Flow, 
Basic  Climatological  Flow 

The  two-dimensional  flow  of  an  inviscid,  incompressible  fluid 
layer,  for  which  the  hydrostatic  approximation  is  made  and  which  is 
subject  to  friction  as  it  moves  over  a  perfectly  spherical  earth  of 
radius  a,  is  governed  by  two  equations  for  the  tangential  components 
of  velocity,  u  and  v: 

(2.1)       ^  u  -  (f  +  H  tan  e)v  +  ^-^^  ^  +   ^^^  =   o   , 

(2o2)    •    ^  V  +  (f  +  I  tan  0)u  +  I  1^  +  Rvq  =  0  , 

where 

_D_  ^  ^       u     ^  +  Z  ^  . 

Dt   "St        a  cos  e   "Sa    a.  W   ' 

and  the  equation  of  continuity  for  the  free  surface  height  of  the 
layer,  h: 

(2.3)        |5  +  ^-^  [^  (hu)  +  ^  (hv  cos  e)]  =  0  ; 

^    -^        dt    a  cos  6    '-dA  ^   ^    d©  ' 

the  system  (2.1)- (2.3)  can  be  represented  as  the  homogeneous  system 

L(W)  =  0  , 

where  W  is  a  vector  with  three  components,  (u, v,h).   The  independent 
variables  A  and  6   represent  longitude  and  latitude  respectively, 
with  corresponding  components  of  velocity,  u  and  v;  g  is  the  gravi- 
tational constant;  f  =  20  sin  0    is  the  Coriolis  parameter,  where  O 


is  the  rate  of  angular  rotation  of  the  earth;  while  R  is  a  constant 
coefficient  so  that  the  friction  force  is  thus  taken  proportional 


r2~.     2 


to  the  square  of  the  velocity,  since  q  =  Ju  +v   represents  the 
speed  of  the  flow. 

In  a  series  of  papers  [5  ;5  ;6  ;  13 ]  the  equations  (2«l)-(2o3) 
are  used  with  R  =  0.   The  quadratic  friction  term  appears  in  the 
Mintz-Arakawa  equations  for  the  motion  of  the  atmospheric  layer 
that  is  next  to  the  surface  of  the  earth  (see  [4  ]).   A  value  for 
R  is  found  by  requiring  that  the  "spin-down"  time  of  the  system 
(2.1)- (2.3)  be  about  10  days.    That  is,  the  variable  W  =  (u,  v,  h) 
is  chosen  at  time  t  =  0  to  represent  a  steady,  purely  zonal  flow 
(see  (2.5))  with  maximum  velocity,  Q;  while  for  the  homogeneous 
difference  equations  that  approximate  the  system  (2.1)-(2.3), 

L^(W^)  =  0  , 

a  value  of  R  is  found  so  that  the  maximum  speed  of  the  solution, 
W.(t),  at  t  =  10  days,  is  reduced  to  .IQ. 

The  order  of  magnitude  of  R  can  be  determined  by  considering 
the  decay  of  the  solution,  u,  of 


This  value  of  the  spin-down  time  was  suggested  by  L,  Gates, 
A.  Kasahara,  C.  Leith,  R.  Somerville,  W.  Washington  and  others  in 
July  1973,  at  an  informal  climatology  meeting  at  the  National 
Center  for  Atmospheric  Research.   Numerical  experiments  on  global 
circulation  models  indicate  that  for  an  atmosphere  which  is 
initially  at  rest,  about  three  weeks  are  required  for  the  normal 
flow  patterns  to  develop  as  a  result  of  the  action  of  the  sun,  etc. 
This  interval  of  time  is  called  spin-up  time.   It  is  natural  to 
call  spin-down  time,  the  time  required  for  the  normal  flow  pattern 
to  dissipate,  in  the  total  absence  of  external  energy  sources. 


^  =  -Ru2  , 
dt 

(2.4) 

u(0)  =  Q   . 

For  Q  =  5  meters/second,  the  value  R  ^  .2(10"-^)  — ^ —  yields 

IIJ c  L/C  X 

u(T)  =  .IQ,  for  T  =  10  daysc   In  fact,  since 


u(T)  = 


Q 
1  +  RQT  ' 


RQT  =  9  is  the  formula  for  determining  R  so  that  u(T)  =  .IQ. 

For  the  frictionless  system  (2.l)-(2.3),  ice.  R  =  0,  a  steady 
zonal  solution  is  given  by: 

u  =  Q  cos  6    , 

(2.5)  V  =  0  , 

2 

h  =  D--|  (aOQ  +  ^)sin^0  . 

The  steady  flow  (2.5)  was  used  to  test  the  difference  scheme  that 
is  "based  on  stereographic  coordinates  (see  Section  3)-   The  fact 
that  the  equations  (2.5)  do  not  depend  on  the  longitude,  makes  the 
zonal  solution  unsuitable  for  fully  testing  finite  difference 
schemes  based  on  "natural"  spherical  coordinates.   Hence,  for  tests 
of  numerical  methods,  it  is  convenient  to  consider  the  flow  equa- 
tions that  result  from  an  "unnatural"  spherical  coordinate  system, 
that  is,  one  fixed  in  the  rotating  "earth",  but  for  which  the  axis 
of  the  coordinate  system  is  perpendicular  to  the  axis  of  rotation 
of  the  earth  (see  [12;13]).   The  new  velocity  variables,  (u,v). 


satisfy  the  system  (2.1)- (2.3)  in  the  new  independent  variables, 
(A,e),  provided  that  the  new  Coriolis  parameter,  f,  defined  by 

(2.6)  f  =  20  cos  A  cos  e    , 

replaces  f  in  (2.1)  and  (2.2).   In  this  case,  the  original  steady 
zonal  flow,  about  the  axis  of  rotation  of  the  earth,  is  given  by: 

u  =  -Q  cos  A  sin  6    , 

(2.7)  V  =  Q  sin  A  , 

2 
h  =  D-i  (aOQ  +  ^)cos^A  cos^S  . 

The  steady  flow  (2.7)  was  used  to  test  the  two  difference  schemes 
that  are  based  on  spherical  coordinates  (see  Sections  4  and  5). 

The  basic  climatological  flow,  W  ,  is  defined  by  using  the 
geostrophic  planar  vortex  flows,  W  ,  studied  by  L.  Bauer  and 
G.  Morikawa  [  2  ] .   That  is,  for  a  rectangular  coordinate  system 
(x,y)  on  a  tangent  plane  that  is  rigidly  attached  at  the  north  pole, 
(O, O),  the  rectangular  velocity  components,  (u  ,v  ),  and  height,  h  , 
of  the  flow,  W  ,  are  given  by 

(2.8)  Up  =  -Yy  , 

(2.9)  -p  =  ^x  ' 

(2.10)  h   =  D  +  -  Y  , 


P 


where,  with  /c  =  — ^ 

/gD 


(2.11)  ^  =  -  W^^^A^^Pi)  • 

X  —  u  i 

Here,  the  variables  (u  ,v  ,ti    )    and  the  stream  function  Y,  are  fionc- 
tions  of  position,  (x,y),  and  time,  t,    by  virtue  of  the  definition, 

(2.12)  ,  p^  S  [X-X.(t)]2+[y_y.(t)]2 

where  (x  (t),y  (t))  is  the  "center"  of  the  polar  low,  while 
(x.  ( t  ),y .  ( t ) )  for  i  =  1,2,3  are  the  "centers"  of  the  semi-permsLnent, 
sub- tropical  highs.   The  constant  7  _>  0  determines  the  strength  of 
the  polar  vortex  while  the  constants  -y-,  =  7p  =  7^  =  -7  "^  0,  deter- 
mine the  strengths  of  the  vortices  that  generate  the  highs.   The 


/  2    2 
distances  j x.   +y.,  are  approximately  the  spherical  distance  between 

the  pole  and  the  sub-tropical  highs  on  the  earth.   The  modified 

Bessel  function  of  the  second  kind  and  order  0,  K  (x),  has  the 

properties : 


K  (x)  +  log  X  -*  0    as    X  -*  0  , 


(2.13) 


K  (x)  r^   l-^  e~^  as   x  — ►  00  . 


The  derivative  of  K  (x)  satisfies: 


(2.14)  K^(x)  =  -K^(x)  , 

where  K-,(x)  is  the  modified  Bessel  function  of  order  1,  for  which: 


8 


K^(x)-i^O     as    x-^0, 

(2.15)  K-^(x)  ^  f^   e"^    as    x  -*  co  , 

K-^(x)  =  -K^(x)  -  i  K^{x)    . 

The  centers  of  the  vortices  satisfy  a  system  of  eight  first 
order  ordinary  differential  equations,  analogous  to  the  equations 
for  the  motion  of  the  logarithmic  vortices  occurring  in  two 
dimensional,  incompressible  fluid  flows  (see  [2  ]).   The  condition 
that  the  vortices  are  in  stationary  equilibrium,  with  the  three 
highs  uniformly  spaced  on  a  circle  of  radius  r,  about  the  polar 
vortex,  is  that 

7       K  (/cryj) 

Under  small  perturbation,  the  system  of  vortices  satisfies  approxi- 
mately a  linearized  system  of  ordinary  differential  equations  with 
constant  coefficients.   The  perturbed  motions  have  several  natural 
periods  of  oscillation  that  depend  on  the  value  of  /cr.   In  fact, 
Bauer  and  Morikawa  find  that  for 

/cr  =  3.75  , 

the  system  has  a  short  period  of  about  1  year  and  a  long  period  of 
about  12  years.   Note  that  for  /cr  =  3  •75^ 


For  /cr  =  3*5,  the  ratio  7  /^  =  .0978...  provides  a  state  of 

stationary  equilibrium;  with  a  short  period  of  about  .8  year  and  a 
long  period  of  about  8  years. 


"Y 

(2.16' )  —  =  .0816...  . 

7 

As  indicated  by  Stewart  [11],  in  his  earlier  study  that  did  not 

D  dArns  s 
include  the  polar  low,  a  mean  sea-level  pressure  of  1.013x10  —^ — p— 

-3  CT  '^"^ 

and  a  constant  density  of  1.22x10   -^-^  ,     correspond  to  a  hydro- 

static  height,  D  ^  8.5km  and  gravity  wave  speed,  /gD  ^  2.87x10  ^^ 


sec 
Hence  for  the  polar  tangent  plane,  rotating  at  the  angular  velocity. 


-  =  >^  ^  1.97  X  lO^km  . 


Therefore 


(2.17)  r  =  r  =  3/75  ^  7,4  X  lO^km  , 

P     Ki 

represents  the  spherical  distance,  r  ,  of  the  sub-tropical  latitude 
circle,  from  the  north  pole.   That  is,  the  sub- tropical  highs  are 
approximately  at  23.5   N  latitude,  in  this  representation. 

Of  course,  the  vortex  approximation  for  the  polar  cyclone  and 
sub-tropical  anti-cyclones  is  singular.   That  is,  from  equations 
(2.10)-(2.13)  it  follows  that 

|h  I  -^  CO    as    (x,y)  ^  (x^,y^)  ; 

furthermore,  the  planar  approximation  uses  a  constant  Coriolis 
parameter,  f.   Therefore,  it  is  a  fact  of  great  interest  that  the 
value  of  r   (with  a  reasonable  choice  of  D)  corresponds  to  the  lati- 
tude  of  the  sub-tropical  highs  as  a  result  of  requiring  the  vortices 
to  have  the  yearly  period  of  the  earth's  revolution  about  the  sun  I 
However,  in  order  to  construct  W  it  is  necessary  to  eliminate  the 


.0 


logarithmic  singularities  and  to  map   the  flows  onto  the  sphere. 

Elimination  of  the  singularities  is  carried  out  first  for  the  planar 

flowo   That  is,  a  planar,  basic  climatological  flow,  denoted  by  W  , 

cp 

is  given  by  equations  (2.8)- (2.11)  with  the  function  K   replaced  by 
the  function  K  ,  defined  as  follows: 
Let 


-J2x 


(2c 18) 


KM  - 


Py(x)  ,  ^  1  ^  1  ^i   ^ 


;(x)  ,  ^i  1  ^   > 


where  x.,  is  some  suitably  chosen  value,  Py(x)  is  a  polynomial  of 
degree  seven,  and  g(x)  is  the  first  term  in  the  asymptotic  expansion 
of  K  (x).   Here,  Py(x)  is  defined  by  requiring: 

P^(o)  =  b,   p^(o)  =  p:^(o)  =  P^"(0)  =  0, 

(2.19) 

P^(x^)  =g(x^),   P^(x^)  =g'(x^),   p:^(x^)  =g"(x^),   P"'(x^)  =  g"'(x^), 


where  b  is  found  by  linear  extrapolation  of  g(x),  from  x  =  x-j^,  i.e, 

(2.20)  b  =  g(x^)  -x^g'  (x^)  . 

In  figure  1,  graphs  of  K  (x)  are  drawn  for  x,  ==  I.762,  2.013,  and 
2.265.   Plots  of  K-l(x),  defined  by 

(2.21)  K^(x)  =  -K^(x)  , 
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appear  In  figure  2,    for  the  corresponding  values  of  x-,  . 

In  other  words,  the  planar,  basic  climatological  flow,  ¥  , 

cp 

is    given  by 

(2.22)  u^p   =  -Oy    , 

(2.23)  ^cp=*x' 

(2.24)  h^^   =   D   +  I  <D    , 


cp 


where 


(2.25)  0    =   -   ^I^I^iKj.p,)     . 

For  a  flow  with  velocity  components  u  and  v,  the  total  relative 
vorticity  in  a  region  R  is  defined  by: 

(2.26)  V(R)  Hjj'  (v^-Uy)dxdy  . 

R 

Integration  by  parts  yields  the  formula, 

(2.27)  V(R)    =(irudx+vdy   , 

BR 

where  the  line  integral  is  evaluated  on  the  boundary,  5R,  of  the 
region  R. 

Consider  the  case  that  the  region  R  is  a  circle  with  center 
at  a  vortex  point,  (x.,y.  ),  and  of  radius  r.  »   Now  the  disteunce 
between  any  two  of  the  vortex  points  is  either  r   or  v^  ^„  (see 
(2.17)).   Hence  for  the  choice 
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it  is  seen  that 


/cr.  =  X,  -  lo762  <  ^p.  , 
X    1  2 


r.  <  -^   (see  (2.17))  . 


It  then  follows  that  for  the  flow  W   generated  by  equations  (2.22)- 

(2.25),  the  derivative  of  the  leading  asymptotic  term  of  K   is  used 

to  calculate  the  total  relative  vorticity  expressed  in  (2.27); 

whereas  for  the  flow  W  generated  by  equations  (2. 8 )- ( 2. 11 ),  the 

derivative  of  the  function  K   is  used  in  the  evaluation  of  the  line 

o 

integral,  (2.27).   It  follows  that,  since  K-,  (x-,  )  and 

K-,(x-,)  =  -g' (x  )  differ  by  less  than  9  percent  for  x-,  =  I.762  ,  the 

total  relative  vorticities,  V(R),  for  the  two  flows  W   and  W 

'       ^     ■"  cp       p 

should  differ  by  about  10  percent  for  the  circle  R  of  radius  r. 
about  (x.,y.).   Hence  the  smoothing  that  produced  W   from  ¥  did 
not  appreciably  change  a  physically  important  feature,  namely  the 
total  relative  vorticity  associated  with  the  cyclonic  and  anti- 
cyclonic  regions  of  the  flow. 

A  description  of  how  the  spherical,  basic  climatological  flow, 

W  ,  is  constructed  from  W   is  given  in  Sections  3-5  for  each  of 
c  cp 

the  three  numerical  schemes  treated  therein.   Suffice  it  to  say  here 
that  if  symmetry  of  the  flow  is  required,  between  the  northern  and 
southern  hemispheres,  then  it  is  necessary  to  smooth  the  planar 
flow,  W  ,  in  the  neighborhood  of  the  "equator".   The  equator  in 
the  polar  tangent  plane  is  taken  at  a  distance  approximately  equal 


K  (x-,  )  differs  from  g(x,  )  by  less  than  6  percent, 
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to  10,000  km  from  the  origin.   Thus  the  formulas  that  represent  a 
"symmetrized",  planar,  climatological  flow  W    are  defined  as 
follows : 


2    2    2     2 
(   u    ,  for  X  +y  =  r   <  d-,  , 

I   cp  "^      —   1  ' 


u 


iq^(r)u^p  +  (1  _^(r))(-^)a,   for   d^  <  r^  <  d^  , 

2    2    2     2 

V    ,  for   y  +  y  =  r   <  d-,  , 

cp  '  "^      —  1  * 

(2-28)   v^^p   ' 

'  ^(r)v^p  +  (l-Q^(r))(|)a,    for   d^  <  r^  <  d^  , 

2    2     2 

h^^  ,  for  X  +y  _^  d-^  , 

2     2     2 
for   d-i  <  r   <  d 
1  —    —  o 


cp  '  J'   _  1 

h    =  I 

scp                        ^  o           o           o 

'  Q^(r)h^p  +  (l-^(r))h  -           ^^     ■    -^     ■    -^ 


In  other  words,  up  to  some  distance  d-,,  from  the  pole,  the 
symmetrized  flow,  W   ,  is  the  same  as  the  planar,  climatological 

o  \_-  L) 

flow,  W   .   Typically  d-,  has  been  chosen  to  approximate  the  spheri- 
cal distance  on  the  earth  from  the  north  pole  to  4  north  latitude, 
while  d   approximates  the  distance  from  pole  to  equatoro   In  this 
band  of  4°  width,  the  planar  flow  is  smoothed  by  interpolation  so 

that  third  derivatives  of  W    are  continuous.   That  is,  %{v)    is  a 

scp  { 

polynomial  of  degree  7  defined  by: 

/■Q^(d^)  =  1  ,    Q^(d^)  =  Q;^(d-^)  =  Q^"(d^)  =  0  , 

(2.29)      \ 

l^(d^)  =  ^(d^)  =  Q^(d,)  =  ^"(d^)  =  0  . 
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It  is  easy  to  verify  that 

0  j<  Qr,(r)  <_  1    for    ^^  —  ^  —  ^o    ' 

The  value  h  is  the  average  value  of  h   evaluated  on  the  circle  of 

cp 

radius  d   about  the  origin;  while,  the  value  u  is  the  average  value 

of  the  tangential  component,  f^^  --^u^  ,  of  the  velocity  (u^p'^cp^ 

evaluated  on  the  circle  of  radius  d   about  the  origino   It  is  easy 

to  verify  that  the  radial  component,  7^^^  t^u^^p,  of  W^^ 

vanishes  at  r  =  d   along  with  all  of  its  derivatives  of  order  less 

than  4.   In  other  words,  W  ^,  is  approximately  a  "zonal"  or 

"circular"  flow  in  the  immediate  vicinity  of  the  "equator",  r  =  d^, 

with  constant  height,  h,  and  with  constant  tangential  component  of 

velocity,  u.   It  is  the  flow  W    that  is  "mapped"  onto  the  northern 
^'  scp 

hemisphere  to  produce  the  spherical,  climatological  flow,  W  ,  by 
symraetrization  across  the  equator. 
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3 .  The  Stereographic  Coordinate  Scheme,  Friction  Coefficient, 
Climatological  Flow 

3a.  Stereographic  Coordinates 

One  of  the  annoying  difficulties,  that  arises  in  the  numeri- 
cal solution  of  hyperbolic  equations  defined  over  a  complete 
sphere,  is  that  any  spatial  coordinate  system  must  have  at  least 
one  singularity.   For  example,  the  spherical  coordinate  system, 
with  lines  of  latitude  and  longitude,  is  singular  at  each  pole. 
The  natural  spatial  mesh  in  spherical  coordinates  uses  uniform 
Intervals  in  latitude  and  in  longitude.   Hence  near  the  poles,  the 
spherical  distance  between  nelghborhlng  mesh  points  is  much  smaller 
than  it  is,  say  near  the  equator.   But,  the  Courant-Frledrlchs- 
Lewy  condition,  that  is  necessary  for  such  a  difference  scheme  to 
be  convergent,  restricts  the  ratio  of  the  time  step  to  the  smallest 
space  Interval.   Thus,  an  Inordinately  small  time  step  must  be 
used,  unless  the  usual  difference  methods  are  altered  in  regions 
where  the  spherical  mesh  distance  Is  small,  i.e.  near  the  poles. 
In  Sections  4  and  ^,    the  use  of  spatial  smoothing  operators  avoids 
this  restriction  of  the  time  step,  at  the  expense  of  additional 
calculations  in  the  neighborhood  of  each  pole.   A  way  to  avoid 
having  such  a  singularity  in  the  coordinate  system,  is  to  cover 
the  sphere  by  two  overlapping  coordinate  patches.   That  is,  a 
spherical  latitude  region  containing  the  northern  hemisphere,  say 
-8   <  0  _<  90  ,  is  mapped  stereographically  onto  a  plane  tangent  at 
the  north  pole,  by  rays  emanating  from  the  south  pole.   An  analo- 
gous stereographic  mapping  is  carried  out  for  a  latitude  region 


17 


containing  the  southern  hemisphere,  say 


i°  >  e  >  -90°  , 


by  projection  from  the  north  pole  onto  a  plane  tangent  at  the  south 
pole.   The  equations  of  motion  are  then  to  be  treated  in  these 
planar  coordinate  systems.   But,  to  avoid  the  natural  singularity 
of  polar  coordinates  in  the  stereographic  planes,  the  spherical 
regions  are  enlarged  so  that  their  stereographic  images  are  squares 
that  are  tangent  to  the  images  of  the  above  described  circles  of 
±8°  latitude,  see  figure  5.   The  stereographic  mapping  is  given  by 

(3.1)        X  =  am  cos  6   cos  A  ,    y  =  am  cos  6    sin  A  , 
where  m,  the  map  factor,  is  defined  by 


(3.2)  m  =  ■= -    -^  <  y  < 


2       _IL  <  e  <  IL 
1  +  sin  B    '  2     -  2  ' 


and  the  radius  of  the  earth  is  a.   The  spherical  velocity  components 
(u,v)  are  defined  by 

(3.3  )  u  =  a  cos  0  A  ,    V  =  a  e  ; 

the  stereographic  velocity  components  {^ijip)    are  defined  by 

(3.4)  ct)  =  i  ,    ^  =  y  , 

where  ( '  )  =  -r-p  (  )  denotes  differentiation  with  respect  to  time. 
By  differentiating  (3.1),  the  relationship  between  the  two  sets  of 
velocity  vectors  is  found: 
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(3.5) 


r  l\ 


ij 


fu^ 


=  -mA 


^^' 


where  m   is    the  map   factor   and  A  is   a  rotation  matrix. 


'       sin   A        cos   A  1 


(3.6) 


A 


-    cos   A        sin   A 


; 


The  equations  of  motion,  (2.l)-(2.3),  become 


(3.7) 


°^        2a  m   ^  ^     ^ 


(3.8) 
and. 
(3o9) 
where 


2a  m  "^ 


-^  H  +  (H(}))^  +  (H^)y  =  0  , 


and  where  the  map  factor  m,  and  Coriolis  parameter  f  =  20  sin  6   are 
expressed  in  terms  of  the  stereographic  coordinates,  and  along  with 


H,  s.  P.,  and  P„  are  defined  by 


(3.10) 


m 


^a^+r^ 


sm 


4 


a 


^4^'   -t^   r2.x2^/  ^ 
4a2  +  r^ 


u  -  h      o2  -  i2  ,  ,2 


m 


l2   ,2 


P^  =   x(({)^  -^^)  +2y4)^  , 
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Since  the  area  of  the  stereographic  image  of  a  spherical  area 

2 
element  is  amplified  by  the  factor  m  ,  the  quantity  H  is  propor- 
tional to  the  mass  per  unit  area  in  the  stereographic  plane, 

3b,  Difference  Scheme 

A  difference  scheme  of  second  order  accuracy  in  both  the  time 
and  spatial  variables  is  now  constructed.  Introduce  a  uniform  mesh 
spacing  in  the  stereographic  square  , 

(3.11)  x^  -  iAx  ,    y  =  jAy 

with 

Ax  =  Ay 

and  let  the  boundaries  be  given  by 

(3.12)  X  =  ±NAx  ,    Y  =  ±NAy  . 
Consider  any  interior  mesh  point 

(3.13)  Pi  =  ^""I'^j^ 

and  its  eight  neighbors  in  counter-clockwise  order, 

P2 "  (^i+i^yj)  '      P3 "  ^^+i'yj+i)  '    P4 "  (^i'^j+i^  ' 
P5  -  (^i-i-yj+i)  '    P6 "  (^-I'^j^  '       P7 "  (^-ryj-i)  ' 

Vq   =    (^i'yj_i)  '      P9  -  ^^i+i'^i-i)  '    ^^®  figure  4  . 


* 


This  is  equivalent  to  using  an  orthogonal  circular  mesh  on  the 
sphere.  That  is,  the  lines  x  =  const.,  y  =  const,  correspond  to 
circles  on  the  sphere. 
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figure  4 

The  points  used  in  advancing  the  solution 
at  the  interior  net  point,  p-,  =    (x.  ,y.) 


Denote  by  L,    ip   ,    R,  m  ,  . . . ,  the  values  of  these  functions  at  time 

t,  at  the  points  p,  ,  ^  Jl  ^  ^  9-   The  scheme  finds  the  values  ({),    , 
^-,    ,  H-,     at  time  t  +At  at  p-,,  by  first  predicting  tentative 

values  at  the  centers  p  ,  p,  ,  p  ,  p ,  of  the  four  squares  that  touch 


P-|_J     i«s. 


(3.14) 


^a.^n   (P1+P2+P5  +P4)    '         Pb   =  ^   (P6+P1+P4+P5) 


Pc   =  ^   (Py+P8+Pi+P6)    ' 


^    (Po  +Pq  +Po+Pt)    • 


For  exajnple,    at   time    t   let 


(3.15) 


K   "^    (Hi+H^+H^+H^) 


and  similarly  define  ^   ,    i^j   ,    etc.   The  predicted  tentative  values 
at  t +At  and  p  ,  i.e.,  ^   ,    tp   ,    Yi     are  found  from  the  Lax-Friedrichs, 
first  order  accurate  scheme.   Let  A  =   At/Ax,  and  set 
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(3.16)  ^^  =  J^-|  [^^[(^2  +  ^3)  -  (^1  +  ^4)]  +iji^^+^i^)  -  (i>2  +  4^i)l 

+  gin^[  (in^H)2  +  (m^H)^  -  (m^H)^  -  (m^H)^,^]  ] 

2a  m  a 

a 

where  s^^  =  <^^  +Tp^; 

(3ol7)   'i'g^  =  Ja"l  Ua^  ^"^2+^3^  -  (^1+^4)]  +?at^^3  )-^4 )  -  (^2  ^^1^^ 
+  gm^[(m^H)^  +  (m^H)^  -  (ni^H)^  -  (in^H)-^]  } 

+  ^M^-faL+-4-  tya^''^a-ia)+2xJ^?J  "iT  Va^"  ' 
2a  m  a 

a 

(3.18)   H^  =  Hg_-|  [[(H^)2  +  (H(t))3  -  (H(t))3_  -  (H^))^] 

+  [i}iTP)^  + {m)i^- {E^y)^- {mp)^])     . 

Analogous  formulas  produce  tentative  values  of  ({)  ,  ^  ,  H  ,  with 
a  =  b,c,d,  at  t  +At  for  the  points  p,,  p  ,  p,.   The  averages  of  the 
four  sets  of  tentative  values,  produce  first  order  accurate,  tenta- 
tive values  at  p.,,  namely 

-Yt+At    1  -^  i 

^1    =  IT  2 ^a   ' 

a=a 

(3.19) 

d 


'Vt+At        1  r 

'}  u 

a 


^1    =  IT  I *-  • 


a=a 


Finally,  the  values  at  t  +At  and  p-,,  i.e.,  (j).    ,    ip-,        '    ^^'^ 
H-,    ,  are  determined  by  a  second  order  accurate,  "corrector"  type 
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formula  (analogous  to  Burstein's  version  of  Richtmyer's  two-step, 
Lax-Wendroff  scheme): 

(3.20) 

+  ^r^^[(^,+^a)-(^c^^b)^^^r^'[^^a^^b)-(*d^^c)^ 

+  gm2[(n.2H^+m2Hj-(m2H^+m2Hj]j 

^  At  v"   f,      ,  7t+AtN  ,    1  r^    fi2       ,2^.7t+Atv2   ,7t+Atx2^ 

+  —  [^1(^1+^1    )+ 2 — t^i^'Pi  -^1  +  (<?!    )  -  (^1    )  ) 

ill  3,  m  -i 


^    o       ^4.    ,      ^A^+At   7t+AtNn         R     /,  .Tt+At   ~t+Atsl 


/~t+At^2         ,Tt+At^2  ^  /7t+At.2 
where    (s^        )      =    (cf,^        )     +  (^1        )    ; 

(3.21) 

+  ^r^'['/',+^a-'/',-^]+^r'''t^a+^-^d-^c] 

+   2x^(^^..^.^fAt  Jt+At)^  __^    (^.s^.V^r^t  '4"^')!      ^ 

(3.22)      H^^^^   =   H^-^    [(H4>)2-(H^)6  +  (Vd  +  Ma)-(Hc4^c-*-«b^b) 

+    (H^)4-(H^)8  +  ^Va  +  Vb)-^Vd  +  Vc)^     • 
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The  values  of  (|),  ^[/,    and  H  at  the  boundary  points  of  the 
northern  stereographic  square  are  found  at  t  +At  by  interpolation 
in  the  values  at  t  +At  for  the  interior  points  of  the  southern 
stereographic  square.   That  is,  the  mapping  of  the  southern  region 
with  rectangular  coordinates  (x,y)  is  defined  by 


(3.23) 


X  =  an  cos  6    cos  A  ,    y  =  an  cos  6    sin  A  , 


where  the  map  factor,  n,  satisfies 


(3.24) 


n 


sm 


1  -  sin  0  ' 
4a2  +7^  ' 


1  <  0  <  IL 


2  - 


2  ' 


-2   -2  ^-2 
r   =  X  +  y 


From  (3.3)^  (3.23 )  and  (3.24),  it  follows  that  the  velocity  com- 
ponents 


(3.25) 

satisfy 

(3.26) 


({)  =  X  ,   f  =  y 


i,  ) 


=  nB 


u 


V 


where  B  is  the  rotation  matrix 


(3.27) 


B  = 


( -    sin  A    cos  A  1 


cos  A    sin  A 


Furthermore,  the  variable 


(3.28) 


i  -  ^ 

n  -  — ^  J 

n 
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is  proportional  to  the  mass  per  unit  area  in  the  southern  stereo- 
graphic  plane. 

The  equations  of  motion,  for  the  southern  stereographic  plane, 
are  identical  to  the  equations  (3o7)-(3.10)  with  (cj),^,  H,  x,  y,m,  f ) 
replaced  by  ((|),^,  H,  x,y,  n, -f ),  but  with  sin  Q   defined  by  (3. 24), 
Hence  the  difference  equations  for  the  advance  of  ~^,    ip ,    H  at 
interior  points  of  the  southern  mesh  (x. ,y. )  are  identical  in  form, 
to  those  for  the  northern  mesh.   But,  each  boundary  point  p  of  the 
northern  square  corresponds  to  some  point  P  on  the  sphere,  whose 
representative  p  lies  well  in  the  interior  of  the  southern  square 
and  vice-versa.   Now  when  a  point  P  on  the  sphere  corresponds  to 
p  =  (x,y)  a  point  of  the  northern  stereographic  square,  then  the 
point  P  corresponds  to  the  point  p  =    (x,y)  in  the  southern  stereo- 
graphic square,  with 

o 
,   ^^,  —   n       n   1  +  sin  Q  4a 

(5.29)  P^HP'    m=T-r^l^  =  -2—2    '    ■ 

Therefore,  if  p  is  a  mesh  point  of  the  northern  plane,  it  is  un- 
likely that  p,  given  by  (3-29).  is  a  mesh  point  of  the  southern 
plane.   In  fact,  for  N  >^  I6,  any  northern  boundary  mesh  point,  p, 
when  plotted  at  its  position  p  in  the  southern  stereographic  plane 
is  surrounded  by  16  interior  mesh  points  of  the  southern  plane,  see 
figure  5.   The  functions  ~^,    f   and  H  are  known  at  the  mesh  points  p^^, 
1  <  k  <  16,  that  surround   p.   The  sixteen-point  Lagrange  inter- 
polation formula  of  the  form 

16    _  _  _ 

(3.30)  "^(P)  =  y~  a^Jp)MP).) 

k=l 
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figure  5 


In  the  southern  stereographic  plane,  the  point  p, 
corresponding  to  a  boundary  point  of  the  northern 
mesh,  is  surrounded  hy  sixteen  interior  points 


is  correct  for  cubic  polynomials «  With  ox  =  0  for  k  =  5^  6,  » .  . ,  l6, 
the  corresponding  four-point  Lagrange  interpolation  polynomial  is 
correct  for  linear  polynomials. 

After  finding  ({)(p),  fiv),    H(^)  by  interpolation,  the  quanti- 
ties ())(p),  ^(p)  are  obtained  from  {3-5)    and  (3.26),  that  is. 


^t^ 


(3o31) 


\ 


^ ; 


-  ^   AB-l 
n 


/  —  N 


i/  ; 


while  H(p)  is  found  from  (3.9)  and  (3.28)  in  the  form 
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2  _ 
(3.32)  H  =  ^  H  , 

m 


where  (3 .29)  is  used  and  it  is  easy  to  verify  that 


AB  ^ 


cos  2A      sin  2A  ; 


sin  2A    -  cos  2A  I 


The  predictor-corrector  type  difference  equations  (3.l6)- 
(3.22)  were  developed  to  solve  the  initial-value  problem  for  the 
homogeneous  system  of  differential  equations  (3o7)-(3«9).   Let  the 
system  of  three  differential  equations  be  denoted  by 

(3.33)  L(W)  =  0  , 

where  W  is  three  component  vector  ((}),(//,  h)  and  L  is  a  vector  with 
three  components,  (L-,,Lp,L,),  that  represent  the  left  sides  of 
equations  (3.7)- (3. 9)  respectively.   Consider  the  corresponding 
non-homogeneous  system 

(3.3^)  L(W)  =  Z  , 

where  Z  has  the  three  components  (E,F,G)  that  are  prescribed  func- 
tions of  (x,y,t).   The  system  (3-3^)  may  be  solved  numerically  by 
using  a  natural  modification  of  the  right-hand  sides  of  formulas 
(3.l6)-(3.18)  and  (3. 20)- (3  .22),  as  follows: 
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In  (3.16),  add  AtE 
in  (3.17),  add  AtF 


in  (3. 18),  add  AtG 


a 


a 


(3.35 


in  (3.20),  addAtE^(t+i^ 
in  (3. 21),  add  AtF-^ 
In  (3.22),  add  AtG-j_ 


(t)    ; 

(t)    ; 

(t)    ; 

(t+^^) 

y 

ft+At^ 

> 

ft+A.^) 

• 

The  solution  of  the  finite  difference  equation  (3*35)  is  denoted 
by  W^. 

3c.  Friction  Coefficient,  Spin-down  Time 

Varying  the  friction  coefficient  R,  will  affect  the  rate  at 
which  a  solution  of  the  initial  value  problem  for  the  homogeneous 
system  (3o7)-(3.9)  (or  equivalently  (3«33))  tends  to  lose  kinetic 
energy.  Numerical  experiments  using  (3 .16)- (3 • 22)  were  performed 
for  initial  data  that  corresponded  to  the  purely  zonal  flow  given 
by  (2.5),  with  a  =  6.37  x  lO^cm;  D  =  8.5  x  10^ cm;  Q  =  50O  ^^, 


sec 


1000 


cm 


•;  R  =  0.5x10  '^  — ,  1.0x10  ^   — .   (For  the  choice  R  =  0, 
sec  cm  cm 

the  purely  zonal  flow  is  an  exact  solution  of  (3 -7 )- (3 '9 ) •  )   The 
region  of  the  northern  stereographic  plane  consisted  of  a  square 
that  circumscribed  the  image  of  the  circle  at  -8   latitude  on  the 
sphere.   The  interval  sizes  were  Ax  =  Ay  =  ^-p  side  of  square 
^  0.916x10^,  At  =  690  sec. 

Figure  6  contains  graphs  of  the  maximum  speed  on  the  sphere, 

I  2    2 

fu  +v  ,  as  a  function  of  the  time,  for  these  four  cases.   For  the 

larger  coefficient  of  friction,  the  maximum  speed  levels  off  after 
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6   Decay  of  the  maximum  velocity  in  zonal  flow 
for  different  friction  coefficients  in  the 
stereographic  coordinates  scheme. 
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cin 
about  5  days,  to  a  value  that  is  about  II5  ;  whereas  for  the 

smaller  coefficient  of  friction,  the  maximum  speed  levels  off 

after  about  6-^    days,  to  a  value  that  is  about  l4o  -.      Both  the 

2    "^  '  sec 

value  of  the  limiting  maximum  speed  and  the  time  at  which  it  is 
attained,  seem  not  to  depend  on  the  initial  speed  but  do  depend  on 
the  friction  coefficient.   The  time  required  to  reach  the  limiting 
maximum  speed  may  be  thought  of  as  being  the  spin-down  time. 

3do  Climatological  Flows 

A  natural  way  to  test  the  accuracy  of  a  long-term  integra- 
tion produced  by  (3-35)  is  to  treat  the  system  [3»3^),    with  the 
non-homogeneous  term  chosen  to  have  the  special  form 

where  W   is  some  prescribed  function.   For  example,  one  such 
choice  that  was  tried  is 


W,  =  (^,,^,,H^)  , 


with 


(bsu    ,    i}j      =   V  ,         mH=h=h    , 

^c     scp  '         '  c     scp  '       c     c     scp 

where  the  symmetrized,  planar  climatological  flow  variables  are 
defined  in  (2.28).   Here  the  centers  of  the  three  anti-cyclonic 
cells  are  initially  placed  120   apart  on  a  circle  of  radius  r   in 
the  stereographic  square  that  corresponds  to  a  latitude  of  about 
30.78°, 


r   =  :x^  +y^  ~  7.24x  10^  . 
p   ^     J   -  I 
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2  2 

Plots  of  the  contour  lines  of  h  ,  H  s  h  /m  ,  H  -  D/m  , 

q  =  ^/ u  +  V  ,  (f)  ,  and  ?//   are  drawn  in  figs.  7,  8a,  8b,  9,    10  and  11 

5 
respectively,  for  the  parameter  values  D  =  8.5x10  cm, 

-y   =  -.08167.,  7.  =  -4x10"^^  cm^/sec  for  i  =  1,2,3;  x.  =  1,762, 

b  =  3.262,  d  =  12.7^x10^,  d-j^  =  Ilo91xl0^,  h  =  8.517x10^ cm, 

u  =  -65.7  where  d  ^  2a  is  approximately  twice  the  radius  of  the 

earth  in  centimeters.   The  solution  W   of  (3*35)  with  initial  data 

W .  =  W  was  computed  for  several  days,  while  the  forcing  function 

Z  =  L(¥  )  was  held  steady,  i.e.  W  was  constant  in  time.   The 

finite  difference  solution  W.  settled  to  a  steady  state  in  about 

2 
30  hours;  a  plot  of  the  steady  contours  of  h  =  m   H. ,  is  shown  m 

fig.  12.   Since  ¥  =  ¥   is  the  solution  of  (3.3'^)^  the  discrepancy 

h.  -  h  is  a  result  of  the  truncation  error  based  on  the  use  of :   a 

A    c 

33x33  spatial  net  point  mesh  for  the  stereographic  square,  the 
corresponding  time  interval  At  =  690  sec,  and  interpolation  to 
find  the  boundary  values.   This  discrepancy  is  analogous  to  the 
initialization  error  that  arises  when  using  real  initial  data  in  a 
global  circulation  model  (see  the  report  IMM-400  by  M.  Ghil  listed 
in  Section  6). 
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Fig.    7 


Contours  of  the  initial  height  field  h^,  plotted  in  the  northern 


stereographic  square, 


The  equator  is  represented  by  the  dashed 
circle.  Maximum  value  of  h^  =  .90x  10  cm  =  M;  minimum  value  of 
h  =  .85xl0^cm  =  m,  contour  interval  =  .554x10  cm.   The  three 

O  Q 

highs  are  centered  at  50  N  latitude  and  are  120   apart  in  longi- 
tude.  The  low  is  centered  at  the  pole. 
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Fig.    8A 


Contours  of  the  initial  height  field  H  ,  plotted  in  the  northern 

2 
stereographic  square  (since  H  = h  /m  ,  the  variation  of  the 

area  factor  m  primarily  determines  the  shape  of  the  contours). 

Maximum  value  of  H  =  .85x  10  =  M,  minimum  value  of 

H  =  642x10-^  =  m,  contour  interval  =  .873x10-^.   The  equator 

is  represented  by  the  dashed  circle. 
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Fig.  8B   Contours  of  the  modified  initial  height  field  H^  -  D/m  ,  with 

/-  o 

D  =    .85x10   ,    plotted   in   the   northern   stereographic    square. 


Maximum  value    of   H    -  D/m' 


.313  X  10' 


M,    minimum  value    of 


H, 


D/] 


m 


102x10^    =  m,    contour   interval  =    .347x10    . 


The 


equator  is  represented  by  the  dashed  circle. 
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Fig.  9  Contours  of  the  initial  wind  speed  q  =  yu  +v  ,  plotted  in  the 


northern  stereographic  square.   Maximum  speed 

,  T  ^      ,  ^4     cm 
=    .116x10     j^       ..,    ^---        -. 33^ 

interval  =    .166x10^     ^^ 


M,  minimum  speed  =  .59x  10"  -^^   =  m,  contour 


sec 


The  equator  is  represented  by  the 


dashed  circle. 
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Fig.  10  Contours  of  the  initial  velocity  component  ^,    plotted  in  the 

northern  stereographic  square.   Maximum  value  of  ^  =  .l6j x  10    =M, 

minimum  value  of  c})  =  -  .IbJ  x  10  =  m,  contour  interval 

=  .571 X  10^^.   The  equator  is  represented  by  the  dashed  circle. 
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Fig.  11  Contours  of  the  initial  velocity  component  f,    plotted  in  the 

northern  stereographic  square.   Maximum  value  of  ^  =  .164x10  =  M, 

minimum  value  of  ^  =  -  . I65 x  10  =  m,  contour  interval 

=  .355xlCr.   The  equator  is  represented  by  the  dashed  circle. 
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Fig.  12   Contours  of  the  height  field  h.,  produced  after  49  hours 

by  solving  the  difference  equations  in  the  stereographic 

coordinates.   Contours  of  the  initial  height  data  h  are 

o 

plotted  in  Fig.  7  at  the  same  interval.   The  truncation 

error  of  the  scheme  produces  the  non-zero  difference  of 


h. 


h 
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4o  A  Mintz-Arakawa  Type  Scheme 
4a.  Difference  Scheme 

The  Mintz-Arakawa  type  scheme  is  based  on  the  conservation 
of  mass  equation  (2.3)  and  the  momentum  form  of  equations  (2.1), 
(2.2);  namely, 

C-l)       Jr  (»>")  *^n^    l^    (»>"')  *^    (hUV  COS  0)] 

+  (f  +-  tan  e)hu  +-i^  45  +  Rhvq  =  0  , 
^    a       '  a  66 

2    2    2 
where  q  =  u  +v  .   A  spatially  staggered,  longitude-latitude  net 

is  used  with  a  simple  predictor-corrector  procedure  to  advance  the 

solution  at  each  time  step.   The  scheme  is  second  order  accurate 

in  the  spatial  increment  and  first  order  accurate  in  the  time 

interval.   This  process  does  not  use  a  cycle  of  five  different 

predictor-corrector  steps  as  is  used  in  the  Mintz-Arakawa  General 

Circulation  Model  [4  ]. 

The  longitude-latitude  mesh  is  based  on  intervals  A?\  =  5  ^ 

A©  =  4°.   The  equator  is  a  latitudinal  mesh  line,  and  the  poles 

are  thus  2°  away  from  their  nearest  latitudinal  mesh  lines.   The 

velocity  components  (u, v)  are  calculated  at  the  intersections  of 

the  mesh  lines.   The  value  of  h  is  computed  at  the  "midpoint"  of 

each  cell  in  the  (A,©)  plane  and  also  at  each  pole.   The  integer 

indices  (i,j)  may  be  used  to  denote  the  points  of  intersection  of 

the  mesh  lines  where  1  <  i  <  I  =  72,  1  <  j  <  J  =  45.   In  this 
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representation  the  points  (1,  j)  and  {7J>,j}    are  identical,  for 

1  <  j  ^  J;  while  the  points  (i,23)  lie  on  the  equator.   The  lengths 

of  the  sides  of  a  mesh  cell  are  given  by 

(4o3)  a  =  a  cos  0a?\  ,         P  =  aAS  , 

in  the  longitudinal  and  latitudinal  directions  respectively.   When 
(i, j)  denotes  the  northeast  corner  of  a  cell,  (u. .,v. .)  denotes 
the  velocity  vector  at  the  northeast  corner  of  that  cell,  while  h.  . 
denotes  the  value  of  the  height  of  the  layer  at  the  "center"  of 
that  cell  in  (a,0)  coordinates  (see  fig.  13 ) • 

It  is  convenient  to  introduce  the  variable 

(4.4)  M.  .   =   aph 

to  represent  the  mass  in  this  cell,  when  both  a  and  h  are  calcula- 
ted at  the  center  of  the  cell.   On  the  other  hand,  the  mass  flux 
per  unit  time  through  the  north  side  of  the  cell  is  denoted  by 

(v.  ,  .  +  V.  .)   (h.  .  +h.  .^-,  ) 

(4.5)  V^.  =  g   "-^-J ^-^ li^^.liJ±J-  , 

where  a  is  evaluated  at  the  latitude  of  the  north  side  of  the  cell; 
while  the  mass  flux  per  unit  time  through  the  east  side  of  the  cell 
is  denoted  by 

(h.  .  +h.^,  .) 

Here,  the  formula  for  the  average  value  of  u.  .,  depends  on 
the  latitude  line  j.   That  is,  if  (u,  ■)^,,^^^„^   were  defined  to  be 
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Fig.  13   The  locations  at  which  the  quantities  with  subscripts 
i,j  are  centered. 
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^  1, J ^average     i, j         2        ' 

then  a  quite  small  time  step  would  be  required  to  ensure  the 
stability  of  the  difference  scheme  (by  the  Courant-Friedrichs-Lewy 
condition).   To  overcome  this  restriction  of  the  time  step,  the 
Mintz-Arakawa  scheme  [ 4  ]  introduces  smoothing  to  compute  the 
fluxes  U^  .,  that  is  the  average  value  of  the  velocity  is  defined 

by 

(N.) 

(4.7)  (u.  .)         =  u.  '^.   , 

^    '^  '  1, J ^average    i, j   ' 

where 


/p,\    u.   .  +  u.   .  -I 
1.  J  2 


i.J     i,J      J   1-1,  J     i,J      1+1,  J         -   -   J 

Here  the  number  of  smoothings  N.  and  the  smoothing  coefficients  /c  . 

are  functions  of  the  latitude  (the  formulas  defining  k-    and  N.  are 

J      J 

given  in  equation  (4.20),  see  [4]).   N.  =  0  for  lines  of  latitude 

J 

near  the  equator.   The  interval  size  At  =  6  minutes  is  now  used. 

Let  the  zero  right-hand  sides  of  the  equations  (4.1),  (4.2) 

and  (2.3)  be  replaced  by  the  functions  F,  G  and  Q  respectively. 

The  difference  scheme  for  predicting  the  value  at  t  +At  of  the  mass 

>•/.  .,  for  2  <  j  <  J  is  found  by  computing  the  mass  flux  across 
1,  J 

the  boundary  of  the  cell  and  the  contribution  from  the  source 
term  Q: 


(4.8)  )^  =  -y       +At 


'    *t  *t  *t  *t  "t    1 

,'"L-l.J-"l,j)*(Vj-l-^i,j'"'"f'«'i,oJ' 
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where  a  and  Q  are  evaluated  at  the  center  of  the  cell  having 
indices  (i,j)  at  its  northeast  cornero   At  the  poles,  h,  and  h   ., 
are  the  values  of  the  height  of  the  layer.   Since  the  surface  area 
-y    of  the  polar  caps  bounded  by  the  lines  of  latitude  corresponding 
to  j  =  1  and  j  =  J  is  given  by 

(4.9)  7  =  ^Tra^  sin^  ^  , 

the  mass  of  the  layer  for  each  of  the  polar  caps  may  be  represented 
by 

(4.10)  ^M^  =  7hi  .   :Vj^i  =  7Vi  • 

The  rate  of  change  of  mass  per  unit  time  in  the  northern  polar  cap 
is  given  by  the  sum  of  the  fluxes  V.  ..  around  the  northern  most 
circle  of  latitude.   Therefore  the  formula  for  predicting  the  value 
of  "M  j^-j_  at  t  +At  is: 

where  a  is  evaluated  on  the  latitude  circle  j   =   Jo   The  corre- 
sponding formula  for  the  southern  polar  cap  is 

(4.12)  M^,^^^   =   J^^^'-Ata  XZ  V*  . 

where  a  is  evaluated  on  the  latitude  circle  j  =  1. 

Formulas  (4.8),  (4.11),  and  (4.12)  can  be  considered  to  be 
of  the  form 

^  t+At  t  ?tV  •,  t 

(4oi3)  ^       =  ir+At[^]    . 
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The  equations  (4.1)  and  (4.2)  (with  right-hand  sides  F  and  G)  are 
next  used  to  produce  predicted  values  of  ^  u  and  ^v  at  t+At,  in 
the  form 

(It. 14) 

That  is,  after  multiplying  (4.1)  and  (4.2)  by  a^,    the  rate  of 
change  — ^--^r — -   and  — ^^^r — -   may  be  determined  when  the  derivatives  -^ 
and  -^   are  replaced  by  -j^  and  -^v-r  ,  respectively  in  the  longitudinal 
and  latitudinal  directions,  as  follows: 

(4.15)    ^%i^  +  ^  (U*u)  +  ^  (V-u)  -aPhv(f  +  ^  ^f  ') 

+   gPh  4^  +  RaPhuq  =  a^F  , 

(/I.l6)         li^  +  ^   (Ifv)   +  ^   (V*v)  +aPhu(f  +  iLlELl) 

rlh 

+  gah  ^  +  RaPhvq  =  aPG  » 


Equation  (4.15)  is  explicitly  centered  at  the  mesh  point 
(i,j)  for  2  <  j  <  J-1,  as  follows: 
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(4.17)  ^  (U"u)  +^  (V-u) 


1 


8      i,j    1+1,  j^^  i,j    1+1,  J    i,J+l    1+1,  j  +  1^ 


i,j    i-l,j'^  i,j    i-l,J    i,J+l    1-1, j  +  1^^ 


[(u.   .+u.  .^.)(V*  .  +V*^,  .  +V*  .^.+V*,^     .,, 
'^  i,j    i,j  +  l'^  i,j    1+1,  J    i,J+l    1+1,  J+1 


-  (u.  .+u.  .  -i)(V*  .+V?,-,  .  +V*  .  -,+V*,.  .  -,)]] 
i,j    i,j-l^^  i,j    1+1, J    i,J-l    1+1, j-1^^^ 


aphv(f  +  ]l_tEL_^)  _  Raphuq  +  a^F 


a 


Xv(f  +  ]L^|£L_£)  _  R;l/uq  +  aPF 


a 


^(^Vi,j^'^i.i,j^'^i,j.i  +  ''Vi.i,j.i) 


u.  .  tan  0  . 
V.  .(f  +   ^^  J  ^ i)  -  Ru,  ,q. 


1.  J 


a 


1,  J  1, J 


+  aPF.  .  , 


where  a  is  evaluated  at  6  .; 


hh  _   gp 


^h> 


^   di    4  '  i,j    1+1, J    i,J  +  l    1+1, j  +  1'  'b± 


average 


where  (^rr-)        is  evaluated  after  smoothing  as  in  (4.7)  through 
Ml  average  ^  \   i  /       t^ 

the  formulas. 
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(Ii) 


(N,) 


"JT  average    ^"ST 


i^n 


(4.18)   (|h)'°'. 


2  L^  1+1,  j+1   i,j+l^   ^  1+1,  J    i,J 


.BhJk-1) 


+  /c 


■^hJk-1)    ,ShJk-l)  ^  ,^hJk-l)' 


1  <  k  <  N.  . 
-   -  J 

The  explicit  centered  formulas  for  equation  (4.l6)  in  the 
range  2  <  j  <  J-1  are  obtained  from  (4.17)  by  interchanging  the 
variables  u  and  v  in  all  the  places  of  (4.17)  for  which  the  corre- 
sponding terms  in  equations  (4.15)  and  (4.l6)  are  so  related,  and 
by  setting 

(4.19)   gah  II  .  g  2^  II  =  ^  UJ^^.  .  -¥,,,,j  +^,,j,,+  ^,,,,j,l) 


The  definition  of  the  smoothing  parameters  /c  .  and  N.  is 

J      J 


based  on  the  formulas 
(4.20) 


^  ^  [§1  ' 


£-1 

a 

8N_. 


for   N.  7^  0  , 


where  the  notation  [x]  represents  the  largest  integer  less  than  or 
equal  to  x,  and  the  latitude  at  which  a  is  defined  depends  on  the 
quantity  that  is  being  smoothed.   That  is,  for  equation  (4.7), 
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since  the  quantity  Ul'  .  is  being  calculated,  a  is  evaluated  at  the 

A  B 
latitude  e  _"   for  2  <  j  <  J;  while  for  equation  (4.18),  since 


hh 

IS  uems  caj.cu±aT:ea,  a  is  evaiuarea  at.  latituae  y 

3 

Furthermore,  on  the  line  j  =  J,  equations  (4.17),  (4ol8),  and 


the  quantity  -;^  is  being  calculated,  a  is  evaluated  at  latitude 


(4.19)  are  used  with  the  variables 

for   1  Jf  i  Jf  I  > 
The  fluxes  U.    -,  are  defined  by  using  (4.6)  with  j    =   J+1  and  with 

(^.22)  (u.  ^^J        =  u   J+1   . 

1,  J+1 'average    i,  j+l 

Here  the  value  u.  '    is  obtained  by  averaging  across  the  pole; 

namely,  with  I  an  even  number, 

u.  ,  -  u   _ 
1 J    .  .  I  -r 

(0)  _       ^  +?^'^ 
^iJ         2        ' 

fk ) 
while  u:  '^  is  defined  by  (4.7)  with  the  parameters  N.  and  /c  . 

1"  3  3 

defined  by  (4.20)  for  a  evaluated  at  0,.   Analogous  formulas  are 

used  for  equations  (4.17),  (4.18),  and  (4.19)  on  the  line  j  =  1. 

By  using  (4.13)  and  (4.l4)  in  this  manner,  predicted  values 

of  i>V.  .  and  then  u.  .  and  v.  .  are  computed  at  each  location  for  the 

time  t  +At.   The  final  corrected  values  are  then  obtained  at  t  +At 

by  using  these  predicted  values  to  approximate  the  derivative  with 

respect  to  time  at  t  +At  that  appears  in  each  of  the  equations: 
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That  is,  the  formulas  that  were  used  to  compute  the  terms  {~tx-)    > 

vf    ,         '^i.^        appearing  in  equations  (4„13)  and  (4„l4),  are  now 

^      ^      ^ 

evaluated  with  the  use  of  the  quantities  ;V  ,  h,  u,  v  at  the  time 

t  +At.   In  [  1  ]  this  use  of  a  backward  time  difference  is  called 
the  Matsuno  time  differencing  scheme. 

4b o  Friction  Coefficient 

The  Mintz-Arakawa  type  difference  scheme  is  more  dissipative 
than  either  the  Lax-Wendroff  type  scheme  described  in  Section  5  or 
the  Kreiss-Oliger  type  scheme  described  in  Section  5-   Hence  for 
the  interval  sizes  AA,  A©  that  were  used,  the  friction  coefficient 
R  for  the  Mintz-Arakawa  type  scheme  is  smaller  by  a  factor  of 
about  10.   In  fig.  l4  and  fig.  15  are  plots  of  the  decay  of  the 
kinetic  energy  of  the  entire  layer  and  of  the  maximum  speed  respec- 
tively, for  the  zonal  flow  of  (2.?),  with  D  =  8.5km,  Q  =  10  -^, 
and  for  each  of  the  values  R  =  5  x  10"^  — ,  8  x  10"-^  — . 

4co  Climatological  Flow 

The  planar  flow,  that  is  described  by  using  the  function 
(I)(x,y)  of  (2.25),  has  the  planar  velocity  components  (iJ.,v)  and  the 
height  h  given  by 
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Decay  of  maximum  speed  in  zonal  flow  for  different  friction 
coefficients  in  the  Mintz-Arakawa  type  scheme. 
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(4.24)  V  =  0^  , 

h  =   D+|$, 

where 

and. 

2    /      ^2  ,  /      n2 


The  height  distribution  could  be  transferred  to  the  sphere  through 
the  mapping  \ 

(4.25)        X  =  a(g-e)cos  A  ,    y  =  a(|-e)sin  A 

and  the  spherical  velocity  components  (u, v)  could  be  defined  by 

fu  =  -|x  sin  A  +  V  cos  A 
-V  =  -\i   cos  A-v  sin  A  . 

In  order  to  simplify  subsequent  calculations,  the  formulas  for  the 
climatological  flow  h,  u,  v  were  defined  on  the  sphere  as  follows 
for  the  northern  hemisphere:  .< 

"5 


(4.27)  h  =  D-^N_^.Kj.p.)  , 


(u,v)  are  determined  from  (4.26)  with  the  quantities  \i   and  v 
replaced  by  ^l  and  v  respectively,  where 
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i=0  p. 

(4.28)  ^ 

1=0  p^ 

Here  the  variables  (x,y)  and  (x.,y.)  are  now  expressed  in  terms  of 
(A,0)  and  (a.,0.)  by  using  (4.25)  and  the  formula  for  p.  is  the 
spherical  distance  between  this  pair  of  points: 

(4.29)  Pi  =  atj). 
where 


CO 


s  (1).  =  sin  6    sin  6.    +   cos  6    cos  0.  cos  (A  -  A.  ) 


Note  that  the  spherical  velocity  field  (u, v)  that  is  hereby 
defined  does  not  correspond  to  the  usual  transformation  of  the 
planar  velocity  field  (ij.,v)  under  the  mapping  (4.25).   That  is, 

f  u  7^  -gr^  A(x,y)  , 


[v   f^^   0(x,y) 

where  (x(t),y(t))  represents  a  particle  path  in  the  planar  flow 
defined  by  (4.24).   Nevertheless,  the  quantities  (u, v)  may  now  be 
considered  to  represent  a  velocity  field  in  the  northern  hemi- 
sphere.  The  essential  cyclonic  or  anti-cyclonic  character  of  a 
planar  low  or  high  is  preserved.   The  flow  is  then  defined  in  the 
southern  hemisphere,  by  requiring  symmetry  across  the  equator. 

Upon  prescribing  [A. ( t ),  0 . ( t )] ,  the  position  of  the  centers 

of  the  subtropical  highs  and  the  polar  low,  the  climatological 

"^     <\^     <\^ 

flow  h,  u,  V  is  defined  as  a  function  of  time. 
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Now,  if  the  differential  operators  appearing  in  equations 
(4.1),  (4»2),  and  (2.3)  are  thought  of  respectively  as 

;r^(u,v,h)  =  0  ,     :- 
:X'^{\x,v,h.)    =   0  , 
,-/   {\i,Y,Y\)    =  0  , 


then  the  forced  system 


X^{u,v,Yi)    =   F  , 


(4.30)  X2^M,Y,lr)    =  G  , 

Z'^tu,  v,h)  =  Q, 

is  solved  by  the  difference  scheme  defined  in  Section  4ao   For  the 
choice, 

F  =     'r-^(u,v,h)  , 

Q  =  x^(u,v,h)  , 

the  system  (4.30)  can  be  represented  as 

(4.31)  /(u,v,h)  =r(u,v,h)  . 


Again,  the  right-hand  sides  are  defined  in  the  southern  hemisphere 
by  symmetry  across  the  equator,  e.g.   v  =  0  on  the  equator.   Hence, 
if  at  time  t  =  0  the  initial  data  are  given  by 


t\j     t\f     *\f  , 


(4.32)  (u,v,h)  =  (u,v,h)  , 
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then  the  equation  (4 ,32)  will  represent  the  solution  of  (4o5l)  for 
all  timeo 

On  the  other  hand,  the  difference  scheme  defined  in  Section 
4a  may  be  thought  of  as  defining  the  functions  (u  ,  v  ,  h  )  that 
satisfy  the  initial  condition  (4,32)  and  the  difference  equation 
system  (on  the  mesh  a) 

(^.33)  ,;r^(u^,v^,h^)  =^(u,v,h)  . 

Hence,  the  errors  u  -u,  v  -v,    and  h  -h  may  be  determined  by  simply 

subtracting  the  known  solution   of  (4.31),  (4„32)  from  the  computed 

values  (u^,v^,h„),  found  in  the  difference  scheme. 
^  A  A   A 

A  plot  of  the  contours  of  the  initial  height  field  h  is  drawn 
in  fig.  l6  for  a  symmetrical  climatological  flow  with  parameter 
values 

D  =  8,429.6m  , 

-y      =  3.26  X  10'''  , 

7i  =  72  =  73  =  -^°0X  10   , 

ic   -  5-886  X  10"^  , 
e^(0)  =^   ,         i  =  1,2,3, 

with  the  function  K  defined  as  in  the  first  graph  of  fig.  1,  for 
X-,  =  l»762o   The  contours  are  plotted  at  an  interval  of  30.28 
meters,  with  the  maximum  value  of  h  s  M  =  8726. 9ni,  attained  at  the 
centers  of  the  highs,  and  the  minimum  value  of  h  s  m  =  8424.2m, 
attained  at  the  poles.   With  this  choice  of  h,  u,  v  held  steady  in 


time,  the  computed  solution  h  ,  u  ,  v  was  found  to  tend  to  a 
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steady  state  within  about  48  hours.   In  fig.  17,  the  contours  of 
h  are  plotted  at  time  t  =  48  hours  at  the  intervals  that  are  used 
in  figo  l6o   Notice  that  there  are  eleven  contour  lines  at  t  =  48 
hours,  in  each  hemisphere,  whereas  there  were  only  nine  at  t  =  0 
hours.   One  of  the  two  new  contour  levels  corresponds  to  the 
maximum  value  of  h  at  t  =  0,  which  appears  since  the  small  trunca- 
tion errors  produce  an  increase  in  the  maximum  value  of  h^  as  time 

A 

increases o   The  additional  closed  contour  lines  that  appear  along 
the  equator  are  evidence  of  slight  elevations,  which  also  arise 
from  the  truncation  error  of  the  scheme. 
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5.  The  Krelss-Ollger  Scheme 

5a.  Difference  Scheme 

A  leap-frog^  explicit  scheme,  that  is  accurate  to  second 
order  in  the  time  increment  and  to  fourth  order  in  the  space  incre- 
ment, had  been  devised  and  used  for  the  system  (2.1)-(2.3)  hy 
Kreiss  and  Oliger  [6  ],  with  R  =  0.   The  variables  (u,v,h)  are 
computed  on  a  uniform  grid  in  (A,0)  coordinates.   Hence,  to  avoid 
the  necessity  of  using  a  small  time  step  At,  Williamson  and 
Browning  [13]  implemented  this  scheme  by  incorporating  the  fast 
Fourier  transform  technique  to  smooth  the  computed  values  of 
(u,v,h)  at  t +At,  on  lines  of  latitude  near  the  poles.   They  were 
kind  enough  to  make  their  program  available.   Subsequently, 
modifications  were  made  here  to  deal  with  the  case  of  forced  motion 
and  friction,  i.e.   R  /  0.   The  details  of  the  scheme  will  now  be 
sketched. 

The  fast  Fourier  transform  algorithm  is  most  efficient  when 
it  uses  a  power  of  2  for  the  number  I  of  net  points  on  a  line  of 

latitude.   Therefore,  the  choices  A0  =  j-jr   corresponding  to  1  =  2  , 

TV  n 

and  AA  =  ^^  were  made.   Let  i//.  .  denote  a  quantity  tp   defined  at  the 

mesh  point  (A., 9.)  at  the  time  t  .   The  generic  partial  derivative 

expressions  i-^  > -^  >  ~^)    appearing  in  equations  (2.1)-(2,3)  are 
replaced  by  the  following  generic  centered  difference  formulas: 

:,,           *n+l   ,n-l 
^^•^  ^  ^  =  ^t^ ^Kl 
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(5.2) 


(5.3) 


hip 


hip 


\^ 

_  4 
3 

h+l,j-^i-l.jl 

1 

3 

2AA 

h^ 

4 
~  3 

>i,j+l-^i,j-l'1 

1 
3 

2A0 

^1+2,-5-^1-2, 


*n+l  /    \ 

Here  if/  represents  the  predicted  value  of  ip   at  time  step  (n+lj. 

Later,  a  formula  for  the  corrected  value  ip  at  time  step  (n+l) 

will  be  given.   The  coefficients  of  the  differentiated  quantities 

in  equations  (2.1)- (2.3)  are  evaluated  at  (A.,0.,t  ).   The  Coriolis 

terms  in  equations  (2.1)  and  (2.2)  are  replaced  respectively  by 


(5.4*) 

where 

(5.5*) 


u 


(f  +-  tan  0)v  , 


*n+l  ,  n-1 
— *    V     +  V 
V   =  o > 


(f  +-  tan  9)'a 
a 


u 


*n+l ^  n-1 
u     +  u 


The  friction  terms  in  equations  (2.1)  and  (2.2)  are  replaced 
respectively  by 


(5.6*) 


Ru^q'^  ,    Rv*q^  , 


with  the  average  values  11^,  v*  defined  in  (5-5*  )•  In  case  the 
right-hand  sides  of  equations  (2.1;-(2.3)  are  non-zero,  say  (E,F,G) 
respectively,  then  these  expressions  are  evaluated  at  (A^,  0.,t^). 

Under  the  difference  prescriptions  (5.I*),  (5.2),  (5.3),  (5.4*)- 

*n+l   *n+l    *n+l 
(5.6*),  the  unknown  predicted  values  u    ,  v    ,  h     appear 

linearly  in  a  set  of  three  simultaneous  equations  with  coefficients 

that  depend  on  quantities  defined  at  the  two  previous  times  t^ 
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and  t^_-|^.   Explicit  formulas  can  thus  be  derived  for  u*'^"^"^,  v*'^"^''", 
h     at  each  net  point  (l,j),  provided  that  special  attention  is 
given  to  the  case  that  (i,j)  is  near  a  pole. 

To  avoid  the  situation  that  (i,j)  is  at  a  pole,  the  values  0. 
were  defined  as  in  the  code  of  Williamson  and  Browning  by 


(5.7)  0j  =  -  g+  (j  --ijAS    for    1  <  J  <  32 


Hence 


ir    A0  ^     ^  TT    A0 


and  (i,j)  does  not  lie  at  a  pole.   It  is  also  apparent  that  (i,j) 

does  not  lie  on  the  equator  either.   When  at  (7\.,0.)    the  value 

1   J 

J  =  1,  2,    ^1,    or  32,  then  one  or  two  points  that  appear  in  the 
difference  expression  5„  will  straddle  the  pole  on  the  same  great 
circle  of  longitude.   In  such  cases,  the  values  of  u  and  v  that  lie 
on  the  other  side  of  the  pole  are  replaced  by  -u  and  -v  respec- 
tively in  (5.3)«   To  prevent  the  separation  of  the  solution  between 
the  even  and  odd  numbered  time  steps,  that  is  characteristic  of 
leap-frog  schemes,  a  corrector  step  is  used  to  define 
(u   ,v   ,h    )  explicitly.   That  is,  equation  (5.1'^)  is  modified 
as 

-.  ,         ,n+l   ,n-l 
^^'^^  ^  =  ^t^ 2At 

the  Coriolis  expressions  {5'^*)    a-nd  (5«5*)  are  modified  as 


(5.4)  (f +^  tan  0)v  ,    (f +-^  tan  e)u  , 

a  a 

where 
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n+1  ,     n-1                         n+1  ,     n-1 
'tri-^                              —        v+v                    —        u+u 
.5.5)^  V   =  5 ,         u  =  ^ 


finally   the    friction   terms    (5 06*)    are   modified  by   using   the 

*n+l 


predicted    speed   q 


(5.6)  Ru^ ^ ,        R?   ^  +q 


■itn+1   ,     n-1                      *n+l  ,     n-1 
■^ — }        R  V  2" 


„    , ,  .     „  ,,     -,     /  n+1  n+1  ,  n+l,  ,  .,  . 

Smoothing  of  the  values  (u   ,v   ,h    j,  on  several  lines 

of  latitude  nearest  the  poles,  is  now  performed  by  the  use  of  the 

fast  Fourier  transform.   This  avoids  the  use  of  a  time  step  At  that 

is  proportional  to  the  arclength  of  the  shortest  mesh  interval, 

i.e.    . 

aAA  cos  (|-^)  ^1  AAA0  . 

The  fast  Fourier  transform  algorithm  determines  the  coeffi- 
cients c,  of  the  finite  Fourier  series,  M-^)?  that  interpolates  the 
2tt    periodic  function  '/'(a)  at  the  points 

Ap  =  p(-^)  ,    0  <  p  <  I  =  64  , 


in  the  complex  form: 

(5o8)  \M  =YZ  %^^^^ 

k=0   ^ 


That  is. 


MAp)  =  ^(Ap)  ,    0  ^  P  1  I  . 


and  the  coefficients  c,  are  uniquely  determined  by  the  sum: 
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-.    I-l       -ikA 

(5.9)  c  =  y  XZ  ^(K^^         ^  '        0  <  k  <  I-l  . 

Observe  that  when  ^(a)  Is  real,  then 

^I_k  "  ~k  '    0  <  k  <  I-l  . 

Hence,  ({)(a)  may  be  written  as  a  real  finite  Fourier  series  in  the 
form 

1/2 

(5.10)  (t>(A)  =  ^    a^  cos  kA  +  bj^  sin  kA  , 


with 


k=0 


a^  =  2  Re  Cj^  ,    b^  =  2  Im  c^^  ,    1  <  k  <  ^  ; 


a^  -  c^  ,    b^  -  0  ;    a^^g  -  ^1/2  '        ^1/2  "  °  ° 

When  the  longitudinal  variable  A  is  used  on  a  line  of  latitude  0, 

277- 
then  the  wavelength  -r-,    of  the  terms  cos  kA  and  sin  kA,  has  the 


arclength  s,  where 


2ir 

s   =  -rr—   a  cos  9    . 
k 


Hence,  the  smoothing  of  i/{'h)  may  be  obtained  by  eliminating  those 
terms  in  (5»10)  for  which  the  arclength  s  is  less  than  — p— •  That 
is,  i/{7\)    is  approximated  by  the  function  l^{7\),    where 

(5.11)  C(^)  =  y^  SL     cos  kA  +  b'  cos  kA  , 

k=0  ^  ^ 


with 


/^x     .  fhir   cos   9      Is 
q  =  q(0)  s  mm  ( — ,  ^) 


=  min  (64  cos  9,    32) 
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Fig.    18      Decay   of  maximum   wind    speed    in    a   zonal    flow,    for 
three   values    of   the    friction   coefficient   R. 
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I 


Here,  if  q  =  _  =  ;52,  then  a^  =  a^,  and  b'  =  b,  for  0  ^  k  j<  q.   But, 
if  q  <  32,  then 


>   _   „  ^,  » 


\   ~    \    '  ^k  "  \    ^°^    0  <  k  <  q-2  ; 


while  for  the  last  two  coefficients 


a"  ,  =  .75a   -,  ,    b'  n  =  .75b   -, 
-       ' -^  q-1  '  q-1     ' -^  q-1 


"q-l  -  "-""-q-l  '  "q-1  '    "-^'^q 

(5.12) 

a  =.25 a    ,     b'=.25a 


For  the  case  I  =  64,  the  values  of  q(0)  were  3,  9^  15j  21,  27  on 

the  five  lines  of  latitude  nearest  to  each  pole.   Hence 

/  n+1  n+1  ,  n+l^         4_,  j  -u   4_     j_ .    j.,„   j.  •      4.  . 
(u   ,v   ,h    )  are  smoothed  by  truncating  the  trigonometric 

interpolation  formula  and  decreasing  the  last  coefficients. 

An  interval  size  At  =  360  sec  was  used  to  insure  stability 

of  the  scheme  for  the  cases  we  treated. 


5b,  Friction  Coefficient,  Spin-down  Time 

The  Kreiss-Oliger  scheme  is  less  dissipative  than  either  of 
the  other  two  schemes.   In  fig,  18,  are  plots  of  the  decay  of  the 
maximum  speed,  that  occur  when  the  initial  state  is  the  zonal  flow 
(2,7)  with  D  =  8.5km,  Q  =  10  -^.   Three  different  values  of  the 

o  C  L- 

friction  coefficient  are  used:   R  =  7.5x10"^  -4r,  5«OxlO"^  -4r> 

■^       cm  cm 

2.5  X  10"^  — c 
-^       cm 

5c .  Climatological  Flow 

Figure  19  shows,  in  stereographic  plane,  the  contour  lines  of 
the  initial  height  h  of  the  layer,  for  the  climatological  flow  with 
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Fig.  19  Contour  lines  of  the  initial  height  h.   The  maximum  level  line  corresponds  to 
h  =  8.66 X  10  cm.   The  minimum  level  line  corresponds  to  h  =  8.46x  10  cm.   The 
contour  interval  is  .02xl0-'cm.   The  letter  M  appears  where  h  ^  8.67x10  cm; 
the  letter  m  appears  where  h  ^  8.44  x  lO-'cm. 
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three  subtropical  highs  at  30   latitude  and  a  polar  low,  as  given 
by  equation  (4.27)  <■   The  mesh  points  used  in  the  scheme  are 
indicated  by  light  dots.   The  value 

-y^   =   -0.4197^  ,    i  =  1,2,3 

was  used  to  ensure  that  the  initial  velocity  at  the  centers  of  the 

highs  is  zeroo   The  other  parameters  in  (4o27)  were  defined  in 

12     2 
c.g.s.  units  as:  -y .    =   -4.5x10   (cm)   for  i  =  1,2,3; 

D  =  8.5 X  10  cm.   The  function  K  was  chosen  as  in  Section  4  (see 

fig.  1  for  X-,  =  1.762),  with  ^rp.  =  3^75^  where  p.  is  the  spherical 

distance  from  the  pole  to  latitude  30  .   The  Coriolis  parameter  in 

(4.27)  is  evaluated  at  the  pole.   The  corresponding  climatological 

velocity  components  u  and  v  are  then  defined  (as  in  Section  4). 

Equations  (2.1)- (2.3)  may  be  written  respectively  in  the  form 


The  forced  system 


'^-^(u,v,h)  =  0  , 


-'r2(u,v,h)  =  0  , 


^-j(u,v,h)  =  0 


V^  /       ,  \      -t^  /'^     '^  O-* 


-^-^(u,  V,  h)  =  Z'^(u,v,  h)  =   E  , 
(5.13)  ^^(u,v,h)  =  <,(u,v,ii)  =  F  , 


^^(u,v,h)  =^^(u,v,h)  3  G  , 


with  initial  condition 

(u,  v,h) 


/  ^V   'V  O./  ' 


u,v,h) 
t=0 
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t=0 


has   the   unique    solution 

(u,  v,h)  =  (u,v,  h) 


r\j     r\,     r\j  ^ 


Note  that  system  (5.I3)  corresponds  to  the  velocity  form  of  the 

equations  of  motion,  while  equations  (4.30)  and  (4.31)  corresponded 

to  the  momentum  form  of  the  equations  of  motion.   The  case  that 

(u, V, h)  is  independent  of  time  was  treated. 

Figure  20  shows  the  contour  lines  of  the  approximate  solution 

h   of  the  system  (5. 13)  at  time  t  =  2  days,  as  obtained  by  the 

predictor-corrector  method  of  Kreiss-Oliger.   The  solution 

(u„,v,,h„)  has  become  steady. 
^  A  A  A 
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Fig.  20  Contour  lines  of  the  height  h.  after  two  days,  produced  by  the  Kreiss- 
Oliger  scheme.   The  contour  levels  are  the  same  as  appear  in  Fig.  19« 
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6.  Summary  and  Publications 

The  three  numerical  schemes  have  been  used  to  compute  non- 
steady  motions  for  periods  of  ten  to  twenty  days.   The  Mintz- 
Arakawa  and  the  Kreiss-Oliger  type  methods  conserve  the  total  mass 
of  the  atmosphere  adequately.   But,  the  stereographic  coordinate 
scheme  ultimately  gains  mass  at  a  small  constant  rate,  after  the 
first  day.   Therefore,  the  stereographic  coordinate  scheme  must  be 
improved  before  it  can  be  used  for  long  term  integrationo   The 
time  step  required  by  this  scheme  is  approximately  twice  as  large 
as  for  the  other  schemes o   A  method  has  been  formulated  for  modi- 
fying a  difference  scheme  so  as  to  force  the  "exact"  conservation 
of  total  mass,  energy,  vorticity,  etc.   This  work  is  not  well 
enough  advanced  to  permit  evaluation. 

Hopefully  long  term  integrations  may  then  be  performed  with 
all  three  schemes,  in  order  (a)  to  decide  upon  the  accuracy  of  the 
average  values  of  such  computed  solutions,  and  (b)  to  determine 
which  type  of  scheme  may  be  the  most  economical  for  application  to 
other  models. 

The  following  list  of  reports,  papers,  and  theses  was  pre- 
pared under  ARPA  sponsorship.   The  thesis  of  M.  Ghil  is  concerned 
with  establishing  the  stability  and  instability  of  the  steady  state 
solutions  of  the  Schneider,  Gal-Chen  energy  balance  climate  model. 

The  thesis  of  A.  Bayliss  studies  the  asymptotic  behavior 
for  t  -^  CO  of  the  solution  of  a  forced  system  of  ordinary 
difference  equations.   The  forcing  function  is  almost  periodic  in 
the  time.   The  difference  equations  result  from  approximating  the 
solution  of  certain  differential  equations. 
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List  of  publications  prepared  under  ARPA  grant 

Bauer,  Lo  and  G.K.  Morikawa,  "Stability  of  rectilinear  geostrophic 
vortices  in  stationary  equilibrium".  Report  no,  IMM-4o6, 
CIMS-NYU,  to  appear  Apr.  1975 » 

Bayliss,  Alvin,  "Exponential  dichotomies  and  almost  periodic  solu- 
tions to  difference  equations",  PhD  thesis  in  preparation, 
NYU,  1975 » 

Drew,  D.A.,  "The  zonally  symmetric  motion  of  the  atmosphere". 
Report  noo  IMM-401,  CIMS-NYU,  Aug.  1973- 

Friedlander,  Susan,  "Spin-down  in  a  rotating  fluid;  part  I", 

Studies  in  App.  Math.,  MIT,  vol.  53,  no.  2,  pp.  III-I36, 

197^0 

Friedlander,  Susan,  "Interaction  of  vortices  in  a  fluid  on  the 
surface  of  a  rotating  sphere",  paper  accepted  by  Tellus, 
1974. 

Ghil,  Michael,  "Numerical  methods  in  fluid  mechanics",  paper  to 
appear  in  Proc.  1973  summer  school  of  Les  Houches,  France, 

Ghil,  Michael,  "On  balance  and  initialization".  Report  no.  IMM-400, 
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Ghil,  Michael,  "A  nonlinear  parabolic  equation  with  applications 
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It  is  intended  that  long  period  integrations  be  performed, 
in  order  (a)  to  determine  the  feasibility  of  obtaining  "accurate 
climatic  forecasts"  by  numerical  methods,  and  (b)  to  evaluate  the 
efficiency  of  the  various  numerical  methods. 

A  brief  summary  is  given  of  the  work  in  progress.   This  is 
followed  by  a  list  of  ARPA  sponsored  CIMS-NYU  publications. 
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